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Free energy calculation of a molecule by removing VDW and Coulomb interactions 
in a transformation and treating the molecule as non interacting systems 
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Free energy calculations in molecular simulations have a variety of applications including deter¬ 
mining the strength of molecular processes such as solvation and binding. It has been recently shown 
that when removing the VDW and Coulomb potential terms of a group of atoms in a molecule by 
performing a transformation, the molecule can be treated as non interacting systems in the free 
energy calculation. This treatment is applicable both when the molecule is in vacuum and in liquid 
and enables a very simple calculation of the free energies associated with the potentials that depend 
on the relative spherical coordinates of these atoms. Here we demonstrate the method in the free 
energy calculation of a Methanethiol molecule in vacuum and in water and compare the results 
to these obtained by MD simulations. The comparison shows agreement between the results and 
faster computation when using the method by factors varying between 5 • 10 3 and 10 12 for the same 
computational resources. 

PACS numbers: 05.10.-a,02.70.Ns,83.10.Rs 


Free energy calculations in molecular simulations are 
used to predict the behavior of the molecules. A variety 
of methods have been introduced both in the context of 
MD (molecular dynamics) and MC (Monte Carlo) M- 
The strength of molecular processes can be estimated by 
comparing the free energy difference associated with the 
process [Bj, 0|. Molecular modeling includes potentials 
that depend on the relative spherical coordinates of the 
atoms such as bond stretching, bond angle and dihedral 
angle potentials and potentials that relate between every 
atom pair in the system such as VDW and Coulomb 
interactions (?]. 

It has been recently shown that when the VDW and 
Coulomb interactions of a group of atoms in a molecule 
are removed in a transformation, the molecule can be 
treated as non interacting systems in the free energy cal¬ 
culation, enabling a simple calculation. This approach 
is applicable both when the molecule is in vacuum and 
in solvent [§]. 

Here we demonstrate the method in the free energy 
calculation of a methanthiol molecule. We then compare 
the results to these obtained by MD simulations in vac¬ 
uum and water environments and compare the running 
time. 


THE DEMONSTRATION 

We calculated free energies of a CH3SH molecule us¬ 
ing the method 0 and compared to the results obtained 
by MD simulations in vacuum/dilute gas and in water. 


We performed the calculations at T = 300K and used 
a realistic force held in which each molecule is individ¬ 
ually parametrized [9] (see Fig. 1) . The force held 
parametrization was according to bond stretching and 
bond angle potentials which are slightly different from 
the standard ones and we performed our calculations 
accordingly. We give the equations for the potentials 
here for completeness: 

V c = (r 2 - t-q) 2 , Vb = ^k e (cos# - cos# 0 ) 2 , 

V d (fajki) = kt/, (1 + cos (n<t> - <t> s )) ■ 


Figure 1. CH3SH molecule with the bond and dihedral angle 
potentials presented, do and kg are on top and bottom and 
are with units of degrees and kj • mol” 1 respectively, r o and 
k c are in units of nm and kj • mol -1 • nm~ 4 respectively. The 
bond stretching terms and their coefficients are ?'oi2 = 0.133, 
k c 12 = 8.87TO 6 , ro23 = 0.183, k c 23 = 5.62TO 6 , ro34 = ^035 = 
ro36 = 0.109, kc 34 = k c 3s = k c 36 = 1-23 • 10 7 . 
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The goal is to demonstrate the free energy calcula¬ 
tions associated with the potentials that depend on the 
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relative spherical coordinates presented in Ref. [§|. This 
calculation is for atoms that do not interact via VDW 
and electrostatic interactions, after that these interac¬ 
tions have been removed in a transformation. Since only 
atoms that are distanced 4 covalent bonds or more in¬ 
teract via VDW and electrostatic interactions, these in¬ 
teractions do not exist in the molecule in the gas state. 
Thus the free energy of the molecule can be decomposed 
into the free energies associated with the bond stretch¬ 
ing terms, the bond angle term, the dihedral angle term 
and the CH3 bond junction. When the molecule is in a 
solvent we removed the VDW and Coulomb interactions 
of the atoms involved, and calculated the free energies 
associated with the bond and dihedral angle terms [ 8 |. 

RESULTS 


—k B T\n2'K + k B T\\\2'Ke ^ k *Io = 

k B T In [e-P k +I 0 (/? k^)] = -0.507982kcal/mol, 

obtained by simply substituting k^ and (3 in this equa¬ 
tion which is based on Eq. (3) in Ref. (g]. This validates 
both the free energy calculation and the partition func¬ 
tion decomposition in Eq. (5). 

We then performed transformations in which the 
bond angle term associated with 9 123 is relaxed. We 
performed the transformation both using the bond an¬ 
gle potential of the force field and by treating the pa¬ 
rameters as if they were given for the standard harmonic 
bond angle term mentioned in Ref. (|j. The results ob¬ 
tained from the MD simulations in vacuum and water 
for the harmonic potential are the following: 


The free energy contributions were calculated by sub¬ 
stituting the force held bond stretching and bond angle 
potentials in the corresponding partition functions in 
Ref. |§] and integrating (see Appendix II) and by using 
the partition functions from Ref. jg]. The results are 
presented below: 

F c i 2 = 5.33kcal/mol, F C 23 = 5kcal/mol, 


AFf, MD vacuum harmonic — 1.39 zfz 0.0023kcal/mol, 

AFfr MD water harmonic = 1.39 zh 0.01kcal/mol. 

These results are in agreement with the calculation ac¬ 
cording to Eq. (2): 

( 2 ) = k B T In2 F^harmonic = 1.38152kcal/mol. 


F c3 4 = F c35 = F c36 = 5.54kcal/mol,F;,i23 = 0.97kcal/mol, 
F d 1234 = —0.58kcal/mol, Fch 3 = 5.39kcal/mol, 


The results obtained from MD simulations in vacuum 
and water for the force held bond angle potential are 
the following: 


where we have used F = — k B TlnZ. These results are AFf, md vacuum cos — —1.388 ± 0.0023kcal/mol, 

valid both for vacuum and solvent as explained in Ref. 


The free energy associated with the bond and dihe¬ 
dral angle terms can be calculated using MD simulation 
by removing terms and subtracting the free energy of 
the corresponding element without the potentials (this 
is simply F = —k B T lnff). The free energy associated 
with the bond stretching terms on the other hand can¬ 
not be calculated by MD simulations since removing 
these terms will result in atoms which are unbound to 
the molecule and therefore impractical calculations. 

To verify the analysis and the calculations we 
performed 21 MD simulations of intermediate sys¬ 
tems which interpolate between the molecule and the 
molecule with the dihedral term associated with tf> 1234 
removed both in vacuum/dilute gas and in solvent. The 
calculated free energy differences associated with these 
transformations were 


AF 6 md water cos = -1.421 ± 0.024kcal/mol. 

These results are in agreement with the calculation ac¬ 
cording to the same potential: 

AFfe cos = —k B T ln2 — Fj, cos = —1.38152kcal/mol. 

We then performed the transformation of the CH3 bond 
junction in which the bond angle terms associated with 
the angles 0234 , # 235 , $ 236 , # 435 , $436 and 0 536 were re¬ 
laxed. The results obtained from MD simulations for 
the harmonic potential are the following: 

AlcH3 MD vacuum harmonic 1=1 9.019 zb 0.007kcal/mol, 

AF(3H3 MD water harmonic = 9.038 zh 0.012kcal/mol. 


AF^md vacuum = —0.5095 ± 0.0023kcal/mol, 

AF^ md water = —0.5095 zh 0.0023kcal/nrol. 
These results can be compared to 


These results are in agreement with result obtained by 
numerically integrating according to Eq. (7) in Ref. [8j 

AFoH 3 Eq. (7) = ~k B T [in (327T 2 ) — In (Zqh 3 harmonic)] — 


AF dEq . (3) = — fee Fin fid — F d 


9.02383 ± 0.0021kcal/mol. 
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We also calculated the free energy difference associ¬ 
ated with the CH3 transformation using the force field 
bond angle potential. The results obtained from MD 
simulations are the following: 

^fcH 3 md vacuum cos = 8.83 A 0.0047kcal/mol, 

AFch3 md water cos = —8.82 ± 0.014kcal/mol. 

These results are in agreement with the result ob¬ 
tained by the numerical integration 

AFch3cos = —8.804 ± 0.00365kcal/mol. 

The calculations of the free energies associated with 
the bond stretching, bond angle and dihedral angle 
terms were immediate ~ 10~ 8 — 10 _7 sec with 8 dig¬ 
its of precision as they are merely substitutions in Eqs. 
(l)-(3) 0 (in the force field we used Eqs. (1) and (2) 
were slightly different). MD simulations in a total run¬ 
ning time of 3.5 • 21 = 73.5 minutes and 1 • 21 = 21 
hours (Six-Core AMD Opteron(tm) processor 2427, 2.19 
Ghz) for vacuum and water environments respectively, 
yielded a precision of 2-3 digits, that may be lower in 
practice due to the spacing between intermediates. 

The running time of the numerical integration in the 
free energy calculation of the CH3 bond junction was 
2.34sec on a 2GHz dual core intel processor. The to¬ 
tal running time of the MD simulations was 3.5 • 21 = 
73.5 minutes and 9 • 21 = 189 hours (Six-Core AMD 
Opteron(tm) processor 2427, 2.19 Ghz) for vacuum and 
water environments respectively. 

DISCUSSION 

We have demonstrated free energy calculation of 
methanethiol by treating the molecule as non interact¬ 


ing systems. The comparison to MD simulation shows 
agreement between the results and faster computations 
when using the method by factors starting from 1800. 
These factors are expected to grow with the molecule 
size as the computation time when using the method 
does not change. 


APPENDIX 

I. SIMULATION PROTOCOL 

We have used the Gromos53a7 force field parameters 
(from ATB server 0) along with spc water model for 
simulations. The cubic box with a minimum distance 
of 1 nm between the solute and the box edge was con¬ 
sidered during the simulations. After minimization, the 
systems were equilibrated (only for solvent simulations) 
in NVT and NPT ensembles for 100 ps. The simulations 
were performed under NPT ensemble, both in vacuum 
and water at each of the 21 equispaced intermediate 
A states including the initial (A = 0) and final states 
(A = 1). The length of simulation was 20 ns in vacuum 
and 1 ns in solvent. However, in the case of bond angle 
transformation of the CH3 group in water, the length 
of the simulation was 20 ns at each intermediate. We 
have computed the free energy change using the BAR 
(Bennett’s acceptance ratio) method in which a series 
of individual free energies is combined into a free energy 
estimate [10]. 


II. ADDITIONAL CALCULATIONS 

We substituted the force field bond stretching and 
bond angle potentials in the corresponding equations 
in Ref. [8(] and integrated. We obtained the following 
expressions 


Z, = 


7re 


8/3k c r 0 


|a (a/4) + Jj (a/4) + a/_i (a/4) + (a + 2) /i (a/4) j , 


where a = f3k c rQ and /„ ( x ) is the modified Bessel function of the first kind, 
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Z nb = / n n e-W cos ^- co <) 2 n sin BidBi dfr 
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